// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_INCOMPLETE_LU_H
#define EIGEN_INCOMPLETE_LU_H

namespace Eigen {

template <typename _Scalar> class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar>>
{
protected:
    typedef SparseSolverBase<IncompleteLU<_Scalar>> Base;
    using Base::m_isInitialized;

    typedef _Scalar Scalar;
    typedef Matrix<Scalar, Dynamic, 1> Vector;
    typedef typename Vector::Index Index;
    typedef SparseMatrix<Scalar, RowMajor> FactorType;

public:
    typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;

    IncompleteLU() {}

    template <typename MatrixType> IncompleteLU(const MatrixType& mat) { compute(mat); }

    Index rows() const { return m_lu.rows(); }
    Index cols() const { return m_lu.cols(); }

    template <typename MatrixType> IncompleteLU& compute(const MatrixType& mat)
    {
        m_lu = mat;
        int size = mat.cols();
        Vector diag(size);
        for (int i = 0; i < size; ++i)
        {
            typename FactorType::InnerIterator k_it(m_lu, i);
            for (; k_it && k_it.index() < i; ++k_it)
            {
                int k = k_it.index();
                k_it.valueRef() /= diag(k);

                typename FactorType::InnerIterator j_it(k_it);
                typename FactorType::InnerIterator kj_it(m_lu, k);
                while (kj_it && kj_it.index() <= k) ++kj_it;
                for (++j_it; j_it;)
                {
                    if (kj_it.index() == j_it.index())
                    {
                        j_it.valueRef() -= k_it.value() * kj_it.value();
                        ++j_it;
                        ++kj_it;
                    }
                    else if (kj_it.index() < j_it.index())
                        ++kj_it;
                    else
                        ++j_it;
                }
            }
            if (k_it && k_it.index() == i)
                diag(i) = k_it.value();
            else
                diag(i) = 1;
        }
        m_isInitialized = true;
        return *this;
    }

    template <typename Rhs, typename Dest> void _solve_impl(const Rhs& b, Dest& x) const
    {
        x = m_lu.template triangularView<UnitLower>().solve(b);
        x = m_lu.template triangularView<Upper>().solve(x);
    }

protected:
    FactorType m_lu;
};

}  // end namespace Eigen

#endif  // EIGEN_INCOMPLETE_LU_H
